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Abstract 

Functions that are smooth but non-periodic on a certain interval pos- 
sess Fourier series that lack uniform convergence and suffer from the Gibbs 
phenomenon. However, they can be represented accurately by a Fourier 
series that is periodic on a larger interval. This is commonly called a 
Fourier extension. When constructed in a particular manner, Fourier ex- 
tensions share many of the same features of a standard Fourier series. In 
particular, one can compute Fourier extensions which converge spectrally 
fast whenever the function is smooth, and exponentially fast if the func- 
tion is analytic, much the same as the Fourier series of a smooth/analytic 
and periodic function. 

With this in mind, the purpose of this paper is to describe, analyze and 
explain the observation that Fourier extensions, much like classical Fourier 
series, also have excellent resolution properties for representing oscillatory 
functions. The resolution power, or required number of degrees of freedom 
per wavelength, depends on a user-controlled parameter and, as we show, 
it varies between 2 and vr. The former value is optimal and is achieved by 
classical Fourier series for periodic functions, for example. The latter value 
is the resolution power of algebraic polynomial approximations. Thus, 
Fourier extensions with an appropriate choice of parameter are eminently 
suitable for problems with moderate to high degrees of oscillation. 

1 Introduction 

In many physical problems, one encounters the phenomenon of oscillation. 
When approximating the solution to such a problem with a given numerical 
method, this naturally leads to the question of resolution power. That is, how 
many degrees of freedom are required in a given scheme to resolve such oscilla- 
tions? Whilst it may be impossible to answer this question in general, impor- 
tant heuristic information about a given approximation scheme can be gained 
by restricting ones interest to certain simple classes of oscillatory functions (e.g. 
complex exponentials for problems in bounded intervals). 

Resolution power represents an a priori measure of the efficiency of a nu- 
merical scheme for a particular class of problems. Schemes with low resolution 
power require more degrees of freedom, and hence increased computational cost, 
before the onset of convergence. Conversely, schemes with high resolution power 
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capture oscillations with fewer degrees of freedom, resulting in decreased com- 
putational expense. 

Consider the case of the unit interval [—1,1] (the primary subject of this 
paper). Here one typically studies the question of resolution via the complex 
exponentials 

f{x) = exp(ia;7ra;), £ M. (1) 

To this end, let 4>n{f), n = 1,2,... be a sequence of approximations of the 
function /(x) = exp(j7rwa;) which converges to / as n — > oo (here n is the 
number of degrees of freedom in the approximation (t>n{f))- For < e < 1, let 
n{e,ijj) be the minimal n such that 

\\f-Mf)\\Li,,^<^- 

We now define the resolution constant r of the approximation scheme {0n} as 

,. ,. n{e,uj) 
r = lim sup lim . 

Note that r need not be well defined for an arbitrary scheme {4>n} (for example, 
if n(e, Lo) were to scale superlinearly in uj). However, for all schemes encountered 
in this paper, this will be the case. 

Loosely speaking, the resolution constant r corresponds to the required num- 
ber of degrees of freedom per wavelength to capture oscillatory behaviour; a 
common concept in the literature on oscillatory problems [6ll20j. In particular, 
we say that a given scheme has high (respectively low) resolution power if it has 
small (large) resolution constant. It is also worth noting that, in many schemes 
of interest, the approximation 4>n{f) is based on a collocation at a particular set 
of n nodes in [—1,1]. In this circumstance, the resolution constant r is equiv- 
alent to the number of points per wavelength required to resolve an oscillatory 
wave (for further details, sec ijl.Sp . 

With little doubt, the approximation of a smooth, periodic function via its 
truncated Fourier series is one of the most effective numerical methods known. 
Fourier series, when computed via the FFT, lead to highly efficient, stable meth- 
ods for the numerical solution of a large range of problems (in particular, PDE's 
with periodic boundary conditions). A simple argument leads to a resolution 
constant of r = 2 in this case (for periodic oscillations), with exponential con- 
vergence occurring once the number of Fourier coefficients exceeds 2ijj. 

However, the situation is altered completely once periodicity is lost. The 
slow pointwise convergence of the Fourier series of a nonpcriodic function, as 
well as the presence of C(l) Gibbs oscillations near the domain boundaries, 
means that nonperiodic oscillations cannot be resolved by such an approxima- 
tion. A standard alternative is to approximate with a sequence of orthogonal 
polynomials (e.g. Chebyshev polynomials). Such approximations possess expo- 
nential convergence, without periodicity, yet the resulting resolution constant 
increases to the value r = tt, making such an approach clearly less than ideal. 

1.1 Fourier extensions 

The purpose of this paper is to discuss an alternative to polynomial approx- 
imation for nonperiodic functions; the so-called Fourier extension. Our main 
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result is to show that Fourier extensions have excellent resolution properties. 
In particular, there is a user-determined parameter T G (l,oo) that allows for 
continuous variation of the resolution constant from the value 2 (in the limit 
T — > 1), the figure corresponding to Fourier series, to tt (T — oo), the value 
obtained by polynomial approximations. Thus, Fourier extensions are highly 
suitable tools for problems with oscillations at moderate to high frequencies. 

Let us now describe Fourier extensions in more detail. As discussed, Fourier 
series are eminently suitable for approximating smooth and periodic oscillatory 
functions. It follows that a potential means to recover a highly accurate approx- 
imation of a nonperiodic function / : [— 1 , 1] — s- M is to seek to extend / to a 
periodic function g defined on a larger domain [— T, T] and compute the Fourier 
series of g. Unfortunately, unless / is itself periodic, no periodic extension g will 
be analytic, and hence only spectral convergence can be expected. Preferably, 
for analytic /, we seek an approximation that converges exponentially fast. 

To remedy this situation, rather than computing an extension g explicitly, 
it was proposed in [7l [11] to directly compute a Fourier representation of / on 
the domain [— T, T] via the so-called Fourier extension problem: 

Problem 1.1. Let G„ he the space of 2T -periodic functions of the form 

n 

g e Gn ■■ g{x) = — + ^ Qffc cos —kx + (ik sin —kx. (2) 

fe=i 

The (continuous) Fourier extension of f to the interval [— T, T] is the solution 
to the optimization problem 

gn ■■= argminll/ - ^11^2 (3) 

Note that there are infinitely many ways in which to define a smooth and 
periodic extension of a function / : [—1, 1] C, of which ^ is but one. We 
say an extension g is a Fourier extension if g G Gn is a finite Fourier series 
on [— T, T]. Note, however, that even within this stipulation, there are still 
infinitely many ways to define g. We refer to the extension gn defined by ^ 
as the continuous Fourier extension of / (in the sense that it minimizes the 
continuous i^-norm). 

As numerically observed in [71 [TT] , the convergence of gn to / is exponen- 
tial, provided / is analytic. When T = 2, this was confirmed by the analysis 
presented in |25j . A pivotal role is played by the map 

y = cos ^x. (4) 

The importance of this map is that it transforms the trigonometric basis func- 
tions that span the space G„ into polynomials iny. In this setting, Problcm ll.il 
reduces simply to the computation of expansions in a suitable basis of nonclassi- 
cal orthogonal polynomials, since the least squares criterion corresponds exactly 
to an orthogonal projection in a particular weighted norm. Well-known results 
on orthogonal polynomials can then be used to establish convergence properties 
of the Fourier extension. We recall this analysis in greater detail in ^ along 
with its generalization to arbitrary T. 

The map ^ demonstrates the close relationship that exists between Fourier 
series and polynomials. Note the similarity to the classical Chebyshev map 
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X = COS 9 that takes Chebyshev polynomials to trigonometric basis functions. 
Compared to Chebyshev expansions, however, it is important to note that the 
roles of polynomials and trigonometric functions are interchanged. The Cheby- 
shev expansion of a given function / is a polynomial approximation, which is 
equivalent to the Fourier series of a related function. The Fourier extension of 
a function on the other hand is a trigonometric approximation, which is equiv- 
alent to the polynomial expansion of a related function. In that sense, Fourier 
extensions and Chebyshev expansions are dual to each other. 

1.2 Key results 

The main result of this paper is that the resolution constant r = r{T) of the 
continuous Fourier extension satisfies 

r(r) < 2Tsin(^) , Te(l,c5o). 

Accordingly, the Fourier extension <;„ of the function ([T]) will begin to converge 
once n exceeds ^r(T)w (recall that Qn involves 2n + 1 degrees of freedom). In 
particular, we find that r(T) - 2T for T w 1 and r{T) - tt as T oo. Note that 
the resolution of the function f{x) = exp(ia;7rx) using classical Fourier series 
on [— T, T] would require a minimum of 2TLtj degrees of freedom (whenever / is 
periodic). Thus, Fourier extensions exhibit comparable performance for T close 
to 1. However, as T increases, / can be resolved more efficiently via its Fourier 
extension than if it had been directly expanded in a Fourier series on [— T, T]. 
In particular, the resolution power is bounded above for all T by tt, which is 
precisely the resolution constant for polynomial approximations. 

Aside from establishing when asymptotic convergence will occur, we also 
show that the convergence in this regime is exponential, and that for all suffi- 
ciently large n the rate is precisely p = E{T), where 

E{T) = cot^ (^) , 

(this result actually holds for all sufficiently analytic functions /, not just ([T])). 
Here, we note that E(l) = 1, implying no exponential convergence for T — I. 
This is of course a consequence of the Gibbs phenomenon of Fourier series on 
[— 1, 1]. Convergence is exponential for all T greater than 1, however, and the 
rate increases for larger T. 

In summary, the main conclusion we draw in this paper is the following. 
Smaller T yields better resolution power of the continuous Fourier extension, at 
a cost of slower, but still exponential, convergence. Conversely, larger T yields 
faster exponential convergence, at the expense of reduced resolution power. For- 
mally, one may also obtain a resolution constant of 2 in the limit n — )■ c» by 
allowing T — )• 1 as n — >■ oo, and towards the end of the paper we shall discuss 
different strategies for doing this, some of which are quite effective in practice. 

Unfortunately, the linear system of equations to be solved when comput- 
ing ^ is severely ill-conditioned. As a result, the best attainable accuracy with 
a continuous Fourier extension is of order -^6, where e is the machine precision 
used. To overcome this, we present a new Fourier extension (referred to as the 
discrete Fourier extension), based on a judicious choice of interpolation nodes, 
which leads to far less severe condition numbers. Numerical examples demon- 
strate much higher attainable accuracies with this approach, typically of order 
e, whilst retaining both the same rates of convergence and resolution constant. 
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Inherent to Fourier extensions is the concept of redundancy: namely, there 
are infinitely many possible extensions of a given function. As we explain, this 
not only leads to the ill-conditioning mentioned above, it also means that the 
result of a numerical computation may differ somewhat from the 'theoretical' 
(i.e. infinite precision) continuous or discrete Fourier extension. Later in the 
paper, we discuss this observation and its consequences in some detail. 

1.3 Relation to previous work 

The question of resolution power of Fourier series and Chebyshev polynomial 
expansions was first developed rigorously by Gottlieb and Orszag [101 p. 35] (see 
also 123]) J who introduced and popularized the concept of points per wavelength. 
Therein the figures of 2 and tt respectively were derived, the latter being gener- 
alized to arbitrary Gegenbauer polynomial expansions in |21| (in 21 provide 
an alternative proof, valid for almost any orthogonal polynomial system (OPS)). 
Resolution was also discussed in detail in [BJ chpt. 2], where, aside from Fourier 
and Chebyshev expansions, the extremely poor performance of finite difference 
schemes was noted. 

One-dimensional Fourier extensions were, arguably, first studied in detail 
in [71 [TT], where they were employed to overcome the Gibbs phenomenon in 
standard Fourier expansions, as well as the Runge phenomenon in equispaced 
polynomial interpolation (see also [5]). Application to surface parametrizations 
was also explored in [TT] . 

Similar ideas (typically referred to as Fourier embeddings) , were previously 
proposed to solve PDE's in complex geometries. These work by embedding the 
domain in a larger bounding box, and computing a Fourier series approximation 
on this domain. See[l[Sl|51[Il[Il[Il[ni[ini[ll[2ilHl[3S]- Such methods 
can be shown to work well, in principle even for arbitrary regions, but they are 
typically prohibitively computationally expensive. 

More recently, a very effective method was developed in [TOl [32] to solve 
time-dependent PDE's in complex geometries. This approach is based on a 
technique for obtaining one-dimensional Fourier extensions, known as the FE- 
Gram method, which is then combined with an alternating direction technique, 
as well as standard FFTs, to solve the PDE efficiently and accurately. This 
approach has also been applied to Navier-Stokes equations [3|. Interestingly, 
it is shown and emphasized in [35] (see also |3]) that this method leads to an 
absence of dispersion errors (or pollution errors) - another beneficial property 
for wave simulations shared with classical Fourier series, and very much related 
to resolution power. 

Having said this, wc note at this point that the FE-Gram approach is quite 
different to the Fourier extensions we consider in this paper. The FE-Gram 
approximation is not exponentially convergent, and in practice, the order of 
convergence is usually limited by the user to ensure stability in the overall 
PDE solver. However, the FE-Gram approximation is designed specifically to 
be computed efficiently, in 0{n\ogn) time, using only function evaluations on 
an equispaced grid. Although there has been some recent progress in the rapid 
computation of exponentially-convergent Fourier extensions of the form ^ from 
equispaced data |30[ I5T] . this is not an issue we shall dwell on in this paper. 
Thus, the main contribution of this paper is approximation-theoretic: wc show 
that one can construct exponentially-convergent Fourier extensions with a reso- 
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lution constant arbitrarily close to optimal. In fJS]we discuss ongoing and future 
work pertaining to these issues. 

This aside, our use of Fourier extensions in this paper to resolve oscillatory 
functions is superficially quite similar to the Kosloff-Tal-Ezer mapping (see [57] , 
[6j chpt 16.9] and references therein, and more recently, [24]) for improving the 
severe time-step restriction inherent in Chcbyshcv spectral methods. Such an 
approach also improves the very much related property of resolution power. In 
this technique, one replaces the standard Chebyshev interpolation nodes with 
a sequence of mapped points, and expands in a nonpolynomial basis defined 
via the particular mapping. Roughly speaking, with Fourier extensions, the 
situation is reversed. Rather than fixing interpolation nodes, one specifies a 
particular basis (i.e. the Fourier basis for a particular T) that gives rapidly con- 
vergent expansions and good resolution properties, and chooses an appropriate 
means for computing the extension (for example, a particular configuration of 
collocation nodes) via the corresponding mapping. 

Finally, we remark in passing that there are a number of commonly used 
alternatives to Fourier extensions. Such methods typically arise from the desire 
to reconstruct a function directly from its Fourier coefficients (or pointwise val- 
ues on an equispaced grid), whether via re-expanding in a sequence Gegenbauer 
polynomials [221 123j , or by smoothing the function by implicitly matching its 
derivatives at the domain boundary [l^. Whilst the latter retains a resolu- 
tion constant of 2, it only yields algebraic convergence of a finite order, and 
suffers from severe ill-conditioning. Conversely, the Gegenbauer reconstruction 
procedure offers exponential convergence, but with a significant deterioration in 
resolution power |21j . 

The outline of the remainder of this paper is as follows. In ^J2]we detail the 
convergence of Fourier extensions for arbitrary T > 1 , and in ^ we address 
numerical computation. In [21 we consider the resolution power of polynomial 
expansions, and in JjSjwc derive the resolution constant for Fourier extensions. 

2 Fourier extensions on arbitrary intervals 

In this section, we recall and generalize the analysis given in [35] for the case 
T = 2. We first show a new result, namely that the continuous Fourier ex- 
tension converges spectrally for all smooth functions / and for all T > 1 fixed. 
Next, a more involved analysis in ^ 2 . 21 demonstrates that the continuous Fourier 
extension does in fact converge exponentially whenever / is analytic. In doing 
so, we repeat some of the reasoning in [55] for the clarity of presentation, as 
well as for establishing notation that is needed later in the paper. 

Before presenting these results, a word about terminology in order to avoid 
confusion. Throughout this paper, we will refer to the following three types 
of convergence of an approximation /„ to a given function /. We say that /„ 
converges algebraically fast to f at rate k if ||/ — /„|| = 0[nr^) as n — > oo. 
Conversely, /„ converges spectrally fast to / if the error ||/ — /n|| decays faster 
than any algebraic power of n~^, and exponentially fast if there exists some 
constant p > 1 such that ||/ — = 0{p~"') for all large n. 
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2.1 Spectral convergence 

Standard approximations based on orthogonal polynomials (or Fourier series 
in the periodic case) converge exponentially fast provided / is analytic, and 
spectrally fast if / is only smooth [161 chpts 2,5]. Whenever / has only fi- 
nite regularity, convergence is algebraic at a rate determined by the degree 
of smoothness: specifically, if / is (fc — l)-times continuously differentiable in 
[—1, 1] and Z*-'^-' exists almost everywhere and is square-integrable (equivalently, 
/ e H'^[—1, 1] - the kth standard Sobolev space of functions defined on [—1,1]), 
then /„ converges algebraically fast at rate k . Our first result regarding Fourier 
extensions illustrates identical convergence in this setting: 

Theorem 2.1. Suppose that f e H'^[-l,l] for some fc e N and that Tq > 1. 
Then, for all n Cz N and T >Tq, 

11/ -. 9„||Lf_,,,, <c..(To)(^)~' 11/11^.,^^^, (5) 

where Ck{Ti^)) > is independent of n, f and T, \\ ■ W^k ^ is the standard norm 

on H^[—l, 1] and gn is the continuous Fourier extension of f on [— T, T] defined 
by ©. 

Proof. Recall that there exists an extension operator £ : H''[—l, 1] H'^{M.) 
with £f\[_ii-^ = / and < c||/||^fc ^ for some positive constant c 

independent of / [T] . 

Let X € C°°(M) be monotonically decreasing and satisfy x(x) = for x > 
To - 1 and x(a;) = 1 for x < 0. We define the bump function B G C°°(M) by 

To<x < -1 
-1 <x <l 
l<x<To 
\x\ > To. 

Note that the function B{x)£f{x) G H^{M.) and has support in [-To,To]. Thus 
its restriction to [— T, T] is periodic for any T > Tq. Since (?„ minimizes the 
norm error over all functions from the set Gn, we have 

\\f-gn\\Ll^^^ < \\f-{B£fULl^^^ < \\B£f-{B£fULi^^^, 

where {B£f)n denotes the nth partial Fourier sum of B{x)£f{x) on [— T, T]. 
Since B{x)£f{x) e H^[-T, T] and is periodic on [— T, T], a well-known estimate 
[M eqn (5.1.10)] gives 

\\B£f-{B£f)n^_^^^^ < {'^y'\\{B£f)^'^h._^^^. 

Moreover, \\{B£ fy'^^h._^_^^ = || (S£:/)(^) < Cfc(ro)||/||^.^^^^^ for some 
constant Ck{To) independent of /. Therefore, the result follows. □ 

2.2 Exponential convergence 
2.2.1 The continuous Fourier extension 

As mentioned in the introduction, the continuous Fourier extension defined 
by ^ can be characterized in terms of certain nonclassical orthogonal poly- 



B{x) = < 



{ xi-x-i) 

1 

X{x- 1) 
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nomial expansions. This was demonstrated for the case T = 2 in [53], but can 
be shown for any T > 1 with relatively minor modifications. 

Our approach is to construct an orthogonal basis for the space Gn of 2T- 
periodic functions. Since the least squares criterion ^ corresponds to an or- 
thogonal projection, it then suffices to expand a given function / in this basis 
in order to find its Fourier extension. 

The cosine and sine functions in Gn are already mutually orthogonal and 
hence can be treated separately. Consider the cosines first, i.e. the set 

Cn {cos —kx}2=Q. 

Trigonometric functions are closely related to polynomials, through an appropri- 
ate cosine-mapping. This can be seen, for example, from the defining property 
of Chebyshev polynomials of the first kind Tk, 

cos kx = Tk (cos x). 

In the same spirit, we define the map 

y = cos ^x, (6) 

and note that, since cosfcx is an algebraic polynomial of degree k in cosx, the 
function cos j^kx is a polynomial in y. From this we conclude that the set of 
cosines C„ is a basis for the space of algebraic polynomials in y of degree n. 

Since the functions in C„ are linearly independent, an orthogonal basis exists. 
Moreover, we may write the basis functions as polynomials in y, say T^{y) — 
T^(cos ^x\ Orthogonality in L^[— 1, 1] implies 



= J J^{cos^x)T^{cos^x)dx 
= 2 ^ Tj(cos ^x)T;^(cos ^x) dx 



= ^ /' Tny)Tf{y)^^dy, 

where, for ease of notation, we have defined the T-depcndcnt constant 

c(T) :-cos^. 

In the latter step above, we applied the substitution (jG]), which maps the interval 
[0, 1] to [c(T), 1] and introduces the Jacobian with the inverse square root. It 
follows that the T^{y) are orthonormal polynomials on [c(T), 1] with respect to 
the weight function 

2T 1 

This weight differs only by a constant factor from the typical weight function 
of Chebyshev polynomials of the first kind Tk{y). However, the interval of or- 
thogonality is different from that of the Chebyshev polynomials, since [c(T), 1] 
is contained within (—1,1] for T > 1, whereas Chebyshev polynomials arc or- 
thogonal over the whole interval [—1, 1]. 
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We now consider the set of sines in C^, 

TT 

Sn {sin-/sa:}'j!^i. 

Analogously, this leads to orthogonal polynomials resembling Chebyshev poly- 
nomials of the second kind Uk{y)- Indeed, from the property 

sin(fc + l)x — Uk{cosx) sin a; 

we find that sin^fcx is also a polynomial in y, but only up to an additional 

factor. This factor is 

sin —X = \J\ - 2/2. 

We therefore consider an orthogonal basis in the form 

Ul {y)V^~^ = C/J (cos ^x) sin ^x. 

Orthogonality in L'^[—l, 1] implies 

/I ^ ^ ^ ^ 

^ f/J(cos -x) sin -xUf {cos -x) sin -xdx 

= 2 / L^J(cos —x) sin —xUf{cos —x) sin —xAx 
_ 2T 

TT 



Ul{y)Uf{y)^l-y^ dy. 

c(T) 



Thus, the Ul{y) are again orthonormal polynomials. The weight function 

W2{y) = — (8) 

TT 

corresponds to that of the Chebyshev polynomials of the second kind, but here 
too the interval of orthogonality [c(r),l] differs from the classical case in the 
same T-dependent manner. 

Since the functions rJ(cos ^x), J7j(cos ^x) sin ^x comprise a basis for the 
space G„, and since the continuous Fourier extension gn is the orthogonal pro- 
jection of / onto G„, we now deduce the following theorem: 

Theorem 2.2. The continuous Fourier extension gn of a Junction f is precisely 

n n — 1 

9n{x) = ^^^k (cos ^x^ + ^ hkUl {cos ^x^ sin ^x, 

k=0 k=0 

where 



J J{x)T^{cos-x)dx, 



afc - / ^ f{xm (cos -x) dx, (9) 

/I ^ ^ 

f{x)Ul (cos -x) sin -x dx, (10) 

and the polynomials T'j^{y) and U"^ {y) are orthonormal on [c(r), 1] with respect 
to the weights wi{y) and W2{y) given by ([7]) and (|H]) respectively. 
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2.2.2 Exponential convergence of the continuous Fourier extension 

The functions being expanded in the theory of the continuous Fourier extension 
above are related to / through the inverse of the map It is necessary to 
distinguish between the even and odd parts of /. 



Since 



a/c = y fix)Tk (cos —x) dx 
= 2 Ux)T^{cos-x)dx 



2T 











lc(T) 





(11) 



we see that the even part of the Fourier extension gn is precisely the orthogonal 
polynomial expansion of the function 

/l(y) = fe COS^l = fe{x). 

A similar reasoning for the odd part of /, based on the coefficients yields 
the second function 

^^^y> ^13^ sin fx' 

with the odd part of gn, when divided by sin ^x, corresponding to the expansion 
of /2 in the polynomials U'^{y). 

Let us now recall some standard theory on orthogonal polynomial (see, for 
example, [SI [351 133 lor an in-depth treatment). The convergence rate of the 
expansion of an analytic function / in a set of orthogonal polynomials is expo- 
nential, and for a wide class of orthogonal polynomials on the interval [—1, 1] 
the precise rate of convergence is determined by the largest Bernstein ellipse 

e{p) ilip-'e-^' + pe'') : 9 e [-tt, tt]}, p > 1, 

within which / is analytic. We shall use the convention p > 1 throughout. Note 
that the ellipse e{p) has foci ±1, and the parameter p coincides with the sum 
of the lengths of its major and minor semiaxes. Suppose now that / is analytic 
within the Bernstein ellipse of radius p — /Omax and not analytic in any ellipse 
with p > Pmax- Then the rate convergence of the expansion of / in a set of 
orthogonal polynomials on [—1, 1] is precisely (pmax)""- 

Returning to the continuous Fourier extension, we now see that its con- 
vergence rate is determined by the nearest singularity of fi{y) or /2(y) to the 
interval [c(T), 1]. Note that, even when / is entire, singularities are introduced 
by the inverse cosine mapping at the points y = ±1. One can verify that the sin- 
gularity at y = 1 is removable for both fi and /2. Thus, the nearest singularity 
lies at y = — 1. 
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To apply the above polynomial theory, it remains to transform the interval 
[c(T), 1] to the standard interval [—1, 1]. This is achieved by the affine map 



™-'(^) = 2^3^-1, .e[c(T),l], (12) 

with inverse m(t) = i(l-c(r))i + 5(l + c(T)) mapping [-1,1] to [c(T),l]. Note 
that maps the point y = — 1 to the point 



c(r) 



(T)-l 



c 



= 1 - 2cosec^ (^) < -1, 



which lies on the negative real axis. Since the Bernstein ellipse e{p) crosses the 
negative real axis in the point — ^(p^^ + p), equating this with u yields the 
maximal value of p. Wc therefore deduce the following result: 

Theorem 2.3. For all sufficiently analytic functions j , the error in approxi- 
mating f by its continuous Fourier extension gn satisfies 

\f{x) - g,,{x)\ < CfE{T)-'\ Vn e N, T > 1, (13) 

uniformly for x G [—1,1], where Cf is a constant depending on f only, and 

E{T) = cot^ (^) . 

This theorem extends Theorem 3.14 of [5S] to the case of general T > 1. 
Note that the estimate holds for all n € N and T > 1. For the sake of 
brevity, we have limited the exposition here to the case of sufficiently analytic 
functions /. Yet, wc make the following comments: 

• The rate of exponential convergence may be slower than E{T) when f{x) 
is not sufficiently analytic as a function of x, i.e. if / has a singularity 
closer than that introduced by the inverse cosine mapping. 

• The convergence may also be faster than E{T) if / is analytic and periodic 
on [— T, T]. In that case, one can verify that the singularity introduced by 
the inverse cosine at j/ = — 1 is also removable. Thus, the convergence rate 
is no longer limited by the singularity introduced by the map between the 
X and y domains. 

• For the case T — 2 we recover the convergence rate E{2) = 3 + 2\/2 found 
in §3.4]. 

The function E{T) is depicted in Fig. [TJ It is monotonically increasing on 
{l,oo) as a function of T, and behaves like 1 + 7r(T - 1) for T « 1 and ifT^ 
for T 3> 1. Thus, larger T leads to more rapid exponential convergence. This 
is confirmed in Fig. [21 where we plot the error in approximating f{x) = e^ 
by /„ for various n and T (for this example we use high precision to avoid 
any numerical effects - see fJS]). Note the close correspondence between the 
observed and predicted convergence rates, as well as the significant increase in 
convergence rate for larger T. Having said this, it transpires that increasing T 
adversely affects the resolution power, a topic we will consider further in fJS] 
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Figure 1: Theoretical convergence rate E{T) of the continuous Fourier extension 
as a function of the extension parameter T. 




(a) Log error (b) Convergence rate 

Figure 2: The quantities e„ = ||/ — 5„|jL~^ (sohd hue) and E[T)^"- (dashed 
hue) for T = 2,4,6,8. The left panel shows the values e„ and E{Ty. The 
right panel gives the scaled values exp(— i loge„) and E{T). 

3 Computation of Fourier extensions 

Although in ij2.2.1l we introduced a new, orthogonal basis for the set G„, this 
basis is never actually used in computations. Instead, we always use either 
complex exponentials or real sines and cosines. The reason for this is primarily 
simplicity (recall that the relevant orthogonal basis relates to nonstandard or- 
thogonal polynomials, which therefore would need to be precomputed), and the 
fact that having a Fourier extension in the form of a standard Fourier scries is 
certainly most convenient in practice. 

With this in mind, let {</>fe}fc"l^ be an enumeration of a 'standard' basis for 
the set Gn (e.g. complex exponentials), and suppose that the continuous Fourier 
extension of a function / has coefficients {afel^"^^ hi this basis. Since (?„ is 
the solution to ([S]) , it is straightforward to show that 

Aa = B, (14) 

where A e c(2"+i)x(2«+i) and B e C^"+^ have entries and 
respectively, and a is the vector of coefficients a^. Thus, given B, we may 
compute the continuous Fourier extension g„ by solving (|14p . 

There are two issues with computing the continuous Fourier extension. First, 
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one requires knowledge (or prior computation) of the integrals {f,(/)j). Second, 
the condition number k{A) ^ E{T)'^'^ Such severe ill-conditioning typically 
limits the maximal achievable accuracy (by this we mean the smallest possible 
error ||/ — ,9n|| that can be attained in finite precision for any n) of the continuous 
Fourier extension to 0{^/e), where e is the machine precision used [25j . 

Ill-conditioning of the exact extension is due to the redundancy of the set 
Goo • This is explained in the next section. Fortunately, its effect can be greatly 
mitigated by computing a so-called discrete Fourier extension instead. This 
extension, which converges at precisely the same rate as the continuous Fourier 
extension, involves only pointwise evaluations of /, and consequently also avoids 
the first issue stated above. This is discussed in ^3.21 



3.1 Ill-conditioning 

The reason for ill-conditioning in the exact Fourier extension is simple: although 
the functions exp(i-^2;) form an orthogonal basis on [— T, T], they only consti- 
tute a frame when restricted to the interval [—1,1]. 

Lemma 3.1. The set 

*:={^}U{-^exp(j^2;)}fegz\{o}, (15) 

is a normalized tight frame for L^[— 1, 1] with frame bound T. 

Proof. Let / e L^[— 1, 1] be given and define / e L'^[—T,T] as the extension of 
/ by zero to [-T,T]. Since 

Ok = f{x)exp{i!^x)dx^—= / /(a;) exp(i^x) dx, (16) 

we find that -^O'k is precisely the fcth Fourier coefficient of / on [—T,T]. By 
Parseval's relation 

oo 

|a,p = r||/f = r||/f, 

k— — oo 

as required. □ 

For a general introduction to the theory of frames, see [T7] . 

Lemma |3. II has several implications. Recall that frames are usually redun- 
dant: any given function / typically has infinitely many representations in a 
particular frame. This implies that there will be approximate linear dependen- 
cies amongst the columns of the Gram matrix A for all large n. This makes A 
ill-conditioned, and in the case of the frame ([15]), leads to the aforementioned 
exponentially large condition numbers. 

Let us consider several particular representations of a smooth function / 
in the frame (fTS)) . Recall that associated to any frame {fk} is the so-called 
canonical dual frame {gk} |17| . A representation of /, the frame decomposition, 
is given in terms of this frame by 



/ = Y{f,9k)fk 



fe=i 
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Since the frame (|15p is tight, it coincides with its canonical dual up to a constant 
factor that is equal to the frame bound, in this case T. Therefore its frame 
decomposition is precisely 

^ J_ exp(i^a;) 

k— — oo 

where the a^'s are given by ((16)) . However, as illustrated in the proof of the 
previous lemma, this is precisely the Fourier series of the discontinuous function 
/, and thus this infinite series converges only slowly and suffers from a Gibbs 
phenomenon. 

On the other hand, as shown in the proof of Theorem 12.11 by extending / 
more smoothly to [— T, T], one obtain representations of / in the frame (fTSj) that 
converge algebraically at arbitrarily fast rates. Thus, we conclude the following: 
different frame expansions of the same function may give rise to wildly different 
approximations. 

Clearly, given its exponential rate of convergence, the continuous Fourier 
extension g„ will not coincide with any of the aforementioned representations 
(recall also that g„ docs not typically converge outside [—1,1] [25], and therefore 
its limit cannot be a frame expansion in general). More importantly, however, it 
is not certain that the solution of the linear system (|14p , when computed in finite 
precision, will actually coincide with the continuous Fourier extension gn itself. 
The reason is straightforward: for large n, A is approximately underdetermined, 
and therefore there will be many approximate solutions to (jl4p corresponding 
to different frame expansions of /. A typical linear solver will usually select 
an approximate solution with bounded coefficient vector a, and there is no 
guarantee that this need coincide with g„. 

This phenomenon has a significant potential consequence: since the numer- 
ically computed extension need not coincide with the continuous Fourier ex- 
tension theoretical results for the behaviour of g„ are not guaranteed to be 
inherited by the numerical solution. Having said this, in numerical experiments, 
one typically witnesses exponential convergence, exactly as predicted by The- 
orem [^31 However, a difference can arise when investigating other properties, 
such as resolution power, as we discuss and explain in ^ 

The above discussion is not intended to be rigorous. A full analysis of the 
behaviour of 'numerical' Fourier extensions is outside the scope of this paper. 
It transpires that this can be done, and a complete analysis will appear in a 
future paper [5]. Instead, in the remainder of this paper, we focus mainly on 
approximation-theoretic properties of theoretical extensions; resolution power, 
in particular. Nonetheless, in ^5.31 we revisit the issue of numerical extensions 
insomuch as it relates to resolution, and give an explanation of the phenomenon 
of differing resolution power mentioned above. 

3.2 A discrete Fourier extension 

The continuous Fourier extension is defined by a least squares criterion, and 
as such it suffers from the drawback of requiring computation of the integrals 
(/, (/ife). To avoid this issue one may instead define a Fourier extension by a 
collocation condition. In other words, if {xi}^2i^ ^ ^'^^ of nodes in [—1, 1] we 
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replace the matrix A and the vector B in by A and B with entries 4)k{xj) 
and f{xj) respectively. We then define Aa = B once more. 

The key question is how to determine good nodes. Recall that as defined 
by © , is essentially a sum of two polynomial approximations in the variable y € 
[c(T), 1]. The polynomial interpolant of an analytic function in y at Chebyshev 
nodes (appropriately scaled to the interval [c(T), 1]) converges exponentially 
fast at the same rate p as its orthogonal polynomial expansion. Therefore, such 
nodes, when mapped to the original domain via x = ^cos~^j/, will ensure 
exponential convergence at rate E(T) of the resulting Fourier extension. 

It is now slightly easier to redefine G„ to be the space of dimension 2n + 2 
spanned by the functions {cos ^kx}^^Q and {sin yfcxj^j^J. The 2n + 2 colloca- 
tion nodes in [—1, 1] therefore take the form {a;i}"^o ^ {~^i}i=0' where 



T -1 
— cos 



i = 0, . . . , n. 

(17) 

Recall that c{T) = cos ^. We refer to such nodes as symmetric mapped Cheby- 
shev nodes, and the corresponding extension g„ as the discrete Fourier extension 
of the function /. 

Besides removing the requirement to compute the integrals (/, ^„), this 
choice of nodes also has the significant effect of improving conditioning, as we 
now explain. We first note the following: 

Lemma 3.2. Let A be the collocation matrix based on the symmetric mapped 
Chebyshev nodes (|17p . and suppose that D is the diagonal matrix of weights 



UJi 



= . Then the matrix Aw = A DA has entries 



iw 



(j)-j{x)(j)k{x)W{x)Ax, fc = l,...,2n + 2. 



where W is the positive, integrable function given by 

2tt coswX 

W[x) 



T yjcos ^x — cos ^ 

Proof. Note that A\y is a block matrix with blocks corresponding to inner prod- 
ucts of sines and cosines with sines and cosines. Consider the upper left block of 
the matrix A^ DA. In the (j, fc)th entry, using the symmetry of the collocation 
points, we have that 

n n 

2^uJi(j)j{xi)(j)kixi) = 2'^ujiTj{yi)Tk{yi), 

i=0 i=0 

where {j/i}"^o = cos(-y;a;i) are Chebyshev nodes on [c(T), 1]. Recall that m{t) is 
the affine map from [—1,1] to [c{T), 1], as defined in H2.2.2[ Since the product 
Tj[y)Tk{y) is a polynomial in y of degree at most 2n, the Gaussian quadrature 
rule associated with the points Xi is exact and it follows that 



f 1 
2y^^ uj^c/) J {xi)(f)k{xt) ^2 / Tj{m{t))Tk{m{t))—==dt 
j-i V 1 — i 



i=0 



1 

T]{y)Tk{y)w{y) dy, 

c{T) 
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(a) Approximation error (b) Condition number 



Figure 3: The approximation errors ||/ — gn\\L^-^ (left panel) and condition 
numbers k{A) (right panel) for the continuous (squares) and discrete (circles) 
Fourier extensions with T ~ 2 and f{x) ~ cos 16a;. 



where w{y) is given by 

w{y) 



1 - <T) _ / ( 2y-l-ciT) 



The first factor is the Jacobian of the mapping to [c{T), 1], the second factor is 
the scahng of the standard Chebyshev weight to that interval. The change of 
variables y = cos now gives 



2 V (a;,)(/)fc {x,) = 2- / 



(x)4>k {x)w{cos ^x) sin ^x dx, 



which is easily found to coincide with {(f>k, 4'j)w ■ The case corresponding to the 
sine functions sin ^kx is identical. □ 

This lemma demonstrates that the normal equations for the discrete Fourier 
extension are the equations of a continuous Fourier extension based on the 
weighted optimization problem. It is well-known that forming normal equa- 
tions leads to worse numerical results Ch.l9], and exactly the same is true 
in this case. Indeed, since W is an integrable weight function, we may expect 
that K,{Aw) ~ E{TY'^ , i.e. exactly as in the unweighted case. Thus collocation 
leads to the significant reduction in condition number, with k{A) « E{T)'^ as 
opposed to E{TY". As a result, one typically observes a much higher accuracy, 
0{e) as opposed to 0{y/e)^ with this approach. 

This improvement is illustrated in Fig.[3]for the example /(x) = cos 16a;. As 
shown in the left panel, the best attainable error with the continuous Fourier 
extension is roughly 10~^, whereas the discrete extension obtains much closer 
to machine epsilon. The right panel indicates the significantly milder growth in 
condition number. 

Such an improvement is by no means unique to this choice of nodes. A 
suitable alternative choice of collocation points follows from using the roots 
of the polynomials T^{y). With an appropriate number of points, the matrix 
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A^DA, with D containing the weights of the associated Gaussian quadrature 
rule on its diagonal, is precisely the matrix of the unweighted continuous Fourier 
extension, rather than the weighted extension that was identified in Lemma 13.21 
The use of these points as Gaussian quadrature was already explored in pS] , 
However, these points depend on T in a non-trivial way and, although they can 
easily be computed, are not available in closed form. As numerical experiments 
indicate that the performance of such point sets is not significantly better than 
the points proposed in this section, we forego a more detailed analysis. 

This aside, we remark that, in some senses, the continuous and discrete 
Fourier extensions are analogous to the orthogonal expansion of a function 
in Chebyshev polynomials and its Chebyshev polynomial intcrpolant (some- 
times known as the discrete Chebyshev expansion). Specifically, the discrete 
versions both arise by essentially replacing an inner product by a quadrature. 
However, the correspondence is not exact, since in the discrete Fourier exten- 
sion the quadrature approximates the weighted inner product on L^[— 1,1]. 
As mentioned, this is done for computational convenience. Nonetheless, the 
approximation-theoretic properties of the discrete and continuous Fourier ex- 
tensions are nearly identical, with the only key difference being that of the 
condition number. 

We have purposefully kept this section on numerical computation of the 
continuous and discrete Fourier extensions short. The main conclusion is that, 
despite quite severe ill-conditioning (even in the discrete case we have exponen- 
tially large condition numbers), one can typically obtain very high accuracies 
(i.e. close to machine precision) with the discrete Fourier extension. Although 
this apparent contradiction can be comprehensively explained, it is beyond the 
scope of this paper, and will be reported in a subsequent work [2]. 

4 Resolution power of polynomial expansions 

Having introduced and analyzed the continuous and discrete Fourier extensions, 
and discussed their numerical implementation, in fj5]we establish their resolution 
power. Before doing so, since similar techniques will be used subsequently, we 
first derive the well-known result for orthogonal polynomial expansions. 

Suppose that f{x) is defined on [—1, 1] and is analytic in a neighbourhood 
thereof. Let /„ € Pn be its n-term expansion in orthogonal polynomials with 
respect to the positive and integrable weight function w{x). If 1, 1] is the 
space of weighted square integrable functions with respect to w, then /„ is given 

by 

/„ := argminll/ -pl|i2 . 
pep„ [-1.11.™ 

In particular, for any p„ G Vn, 

Hence, we note that to study of the resolution power of any orthogonal poly- 
nomial expansion it suffices to consider only one particular example. Without 
loss of generality, we now focus on expansions in Chebyshev polynomials (i.e. 
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w{x) (1 — .T^) 2)^ in which case 

n 00 

Pnix) =^ 'akTk{x), f{x)^^'akTk{x), (18) 
fc=o fc=o 

where ^ 

Ofc = / /(cos6')cosfc6id6', (19) 



and ' indicates that the first term of the sum should be halved. Since the infinite 
sum (jl8|) converges uniformly, we have the estimate 

||/-Pn||Lp,_„ < El^'^l- 

Therefore it suflices to examine the nature of the coefficients a„ for the function 

f{x) = exp(i7rwa;). (20) 
We use the following standard estimate, given in [35l p. 175]: 

Wn < 21 

Here p corresponds to any Bernstein ellipse e{p) in which / is analytic and Mp 
is the maximum of 1/(2) | along that ellipse. 

For a given uj and n, we consider the minimum of all bounds of the form (|2ip . 
Since / is a complex exponential, / reaches a maximum at the point on the 
ellipse e{p) with the smallest (negative) imaginary part. This corresponds to 
6 — —it/2 and thus we have 

Mp = exp (-KLo^^-—^ . (22) 



The bound (1211) now becomes 



2p 



pU pTl 

For fixed uj and ri, we find the minimal value of this bound by differentiating (|23p 
with respect to p. Equating the derivative to zero 

-^B(c^,n,p) = 0, (24) 
dp 

yields two solutions 

n ± ^V? — TT^W^ 

■ 25 

TTo; 

Consider the case n < ttlo. Both solutions of ([24)1 are complex- valued. Note 
that 

B{Lj,n, 1) = 2, 

and 

d 

—B[uj, n, 1) = 2710; — 2n. 
dp 
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It follows that B is strictly increasing as a function of p. For n < ncu, the best 
possible bound for a„ of the form (PT|) is 2. 

Consider next the case n > ttuj. Since it is easily seen that the roots are 
inverses of each other, we restrict our attention to the one greater than 1, i.e., 



Pn 



Since B{uj,n,p) initially decays at p = 1, pmin is the unique minimum of the 
bound for p > 1. Thus, 

B{uj,n,p,nin) < B{uj,n,p), pe[l,oo). 

In particular, this implies exponential decay of the next coefficients: 

1 2 

|a„+fc| < B{uj,n, p,nin)— — < — — ■ 

Pmin Pmin 

Finally, consider the value n = ttcj. Then p = 1 is a global minimum of the 
bound (l23l). since 



B{ijj, noj, 1) = -y^B{u}, TTUJ, 1) = 0, 



-i3(c.,™,l) = - 
and 

-—rB(uj, TTUJ, 1) = 2tTUJ > 0. 

dp'^ 

In conclusion, we have seen that for an oscillatory function of the form (PO)) . 
we need n = ttuj degrees of freedom before exponential decay of the coefficients 
an sets in. This corresponds to tt degrees of freedom per wavelength - the 
well-known result on the resolution power of polynomials. 

Note that this value was previously derived in |21| for orthogonal polynomial 
expansions corresponding to the Gegenbauer weight w (a;) = (1— a;^)'^"^^ A > — ^ 
(see also [20l p. 35] for A = 0). This result was based on explicit expressions 
for the Gegenbauer polynomial coefficients of e'^*. The previous arguments 
generalize this result to arbitrary weight functions w{x). In fact, we could have 
also used the explicit formula for the Chebyshev polynomial expansion of e'^* 
to derive estimates similar to those given above (and possibly more accurate). 
However, we shall use similar techniques to those presented here in the next 
section, where such explicit expressions are not available. 



5 Resolution power of Fourier extensions 

We now consider the resolution power of the continuous and discrete Fourier 
extensions. As commented in 331 the continuous/discrete Fourier extension <;„ 
need not be realized in a finite precision numerical computation. Hence, we 
divide this section between theoretical estimates for the resolution power of the 
continuous/discrete extension, and its numerical realization. 
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5.1 Resolution power of the continuous Fourier extension 

A naive estimate for the resolution constant r(T) follows immediately from the 
bomrd ^ in Theorem 12. II Indeed, for f{x) ~ exp(ia;7rx) we have 

11/11^.^^^^^ = /c = i,2,..., 

and therefore r(T) satisfies r(T) < 2T, with spectral convergence occurring once 
n exceeds ujT. 

On first viewing, this estimate seems plausible. For example, consider the 
special situation where the frequency of oscillation uj is an integer multiple 
of T~^. Then the function f{x) = exp(i^mx) is precisely the mth complex 
exponential in the Fourier basis on [— T, T]. Given that the continuous Fourier 
extension gn of / is error minimizing amongst all functions in G„ , we can expect 
/ to be recovered perfectly (i.e. gn = /) by its continuous Fourier extension 
whenever n > m = loT. Thus, for oscillations at frequencies a; = ^, m G Z, the 
estimate r{T) < 2T appears correct. 

However, empirical results indicate that such an estimate is only accurate 
for small T: for large T it transpires that r(T) ^ tt. We shall prove this result 
subsequently. Before doing so, however, let us note the following unexpected 
conclusion: when T is large, we can actually resolve the oscillatory exponentials 
exp(iyTOa;) accurately on [—1, 1] using Fourier extensions comprised of relatively 
non-oscillatory exponentials exp(iyna;), n <^ m. 

To obtain an accurate estimate for r(T), we need to argue along the same 
lines as ^^2.21 and exploit the close relationship between Fourier extensions and 
certain orthogonal polynomial expansions. Recall that the theory of ^2.21 treats 
even and odd cases separately. Let us first assume that / is even, so that 

f{x) ~ cosojTra; 

(we consider the odd case later). Upon applying the transformation x — 
^ cos""'^ y we obtain 

/i(y) ^cosL^Tcos-iy, 2/G[c(r),l]. 



The Fourier extension of f{x) is precisely the expansion of fi{y) in the orthog- 



onal polynomials T^{y). If we now map the domain [c{T), 1] to [—1, 1] via 



c{T) + l-2y u , , irrw 

^= c(r)-l ' V-l^{l-c{T)) + -{l + c{T)l 

then this equates to an orthogonal polynomial expansion of the function 



/2(m) = cos 



c.rcos-1 [^{l-c{T)) + \{l + c{T))^ 



(26) 



Thus, as in 21 to determine the resolution power of gn, it suffices to consider 
the expansion of /2(u) in Chebyshcv polynomials on [—1, 1]. 

In view of the bound (j2ip . we now seek the maximum value of /2(u) along 
the Bernstein ellipse e{p). We first require the following lemma: 

Lemma 5.1. For p < E{T) and sufficiently large uj, the maximum value of 
1/2(^)1 on ttie Bernstein ellipse e{p) occurs at 9 = 0. 
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Proof. Let z ~ X + iy. Then 

I coswzp = ^ (cos2a;a; + cosh2ajy) . 

Hence, for sufficiently large uj, the maximal value of | cos ijJz\ on some curve C 
in the complex plane occurs at the point z € C where is maximized. 

When u = i(p-ie-'* + pe'*), we may write /2(m) = cosujT[X{0) + iY{9)], 
where X{d) and Y{9) are defined by 

u 1 

cos [X{0) + iY{e)] - -(1 - c(r)) + -(1 + c(T)) = A(^) + ii?(0), 
and, letting p = for 77 > (recall that p > 1), 

A(e') = i(l-c(r))cosh?7Cos6l+i(l + c(r)), B{0) = ^(1 - c(r)) sinh7]sin6'. 

Expanding the cosine, and equating real and imaginary parts, wc find that 

cosX(6') coshy(6l) = A{e), amX{e) sinhy(6l) = B{9). 

We seek to maximize Note trivially that Y{9) cannot vanish identically 

for all 9, and thus it suffices to consider only those 9 for which Y{9) > 0. Let 
Z{9) = cosh^ Y{9). Then Z{9) is defined by 

A'i9) B\9) 
Z{9) Z{9) - 1 

Rearranging, 

Z'^{9) - (1 + A^{9) + B^{9))Z{9) + A'^{9) = 0. 

To complete the proof, it suffices to show that Z attains its maximum value at 
9 = 0. Suppose not. Then Z'{9o) = for some 6*0 7^ (with Z{9o) > 1), and, 
after differentiating the above expression, we obtain 

{A\9o)A{9o) + B\9o)Bi9o)) Z{9^) = J^(9^)A(9^), 

i.e. 

sin 6*0 ((1 - c(T)) cos 6*0 +(1 + c(T)) cosh ry) ^(6*0) 

= sin 6*0 cosh?? (1 + c(T) - (c(T) - 1) cosh cos 6*0) . 

Let us assume first that 9q ^ 0, ±7r. Hence 

((l-c(T)) COS6I0 +(l + c(T)) cosh 77)^(6*0) 

= cosh?? (1 + c(T) - (c(T) - I)cosh7?cos6'o) . (27) 

Suppose the left-hand side vanishes of (P7|) vanishes, i.e. 

cos 00 = ^77^;^ cosh?/. (28) 
c(T) - 1 
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Since sin(0o) 7^ the right-hand side of (P7|) must also vanish, and therefore 

(29) 



cos 6*0 = ^TT^^SiZlgechTy. 



c(r) - 1 

Equations ([28| and ((29|) cannot hold simultaneously (since 77 > 0), and therefore 
the left-hand side of (|27p cannot vanish. We deduce that 



Zi0o) = 



cosh?7(l + c(r) - (c(T) - I)cosh77cos6lo) 



(1 - c(r)) cos6'o + (1 + c(r)) cosh?; 
We now substitute this into the quadratic equation for Z(6) to give 

(cos 6I0 - cosh r;)^ ((c(T) - 1) cosh r; cos - {1 + c{T))) = 0. 
Since 77 > 0, we must have that 

1 + c{T) 



cos 9q 



-secliTy. 



(c(T) - 1) 

However, substituting this into (|27| . gives Z{9o) = 0, a contradiction. Therefore 
6*0 = 0, ±7r. It is easy to see that 6*0 = ±7r leads to a smaller value of Z{do) than 
6*0 = 0. Hence the result follows. □ 

Using this lemma, we deduce that the 7T.th coefficient a„ of the Chcbyshev 
expansion of /2(m) is bounded by 



B{uj,n,p,T) = — cos 



ujT COS 



'"Up+-) (l-c(T)) + i(l + c(T)) 



■ cos 



p" 



Tcosh-i Q l^p + ^) (1 - c{T)) + 1(1 + c{T)) 



We proceeded in g] by computing the roots of the partial derivative of the 
bound with respect to p. The same approach applies here, but unfortunately it 
does not lead to explicit expressions. However, a simple modification makes the 
bound more manageable. We write 

cos(x) = i(e*"+e-*"). 

Since the argument of the cosine for p > 1 is purely imaginary, with positive 
imaginary part, the second exponential dominates the first. Hence, for large oj, 
we may approximate B(w, n, p, T) by 



B{uj,n,p,T) 
1 



pn 



c.rcosh'1 ( 1 (^p + (1 _ c{T)) + 1(1 + c{T)) 



(30) 



One can then explicitly find roots of the partial derivative of B with respect to 
p. We have 

p*in) 



\c{T) + 3) + io^T^{c{T) - 1) ± 2n^w2T2(c(T)2 - 1) + 2n^{c{T) + 1) 



w2r2(c(T)-l) + n2(l-c(r)) 



(31) 
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These two roots are again each other's inverse. One quickly verifies that the 
square root is real and positive if and only if 

n>ic^rV2-2c(T)-ic^r(r), (32) 

with 

r{T) :=Tv/2-2c(r) = 2rsin(^). (33) 

If the condition (|32|) is satisfied, selecting the + sign in ((3T|) yields the root that 
is greater than 1. 

A little care is necessary when applying the bound (15(11) . Recall that /2(w) 
is only analytic in Bernstein ellipses e{p) with p < E{T). It may be the case 
that p*{n) > E{T) and therefore we cannot use this bound directly. However, 
explicit computation yields the condition 

2 1 

E(T) > p*(n) ^ n> ujT = =ujT. (34) 



■.2T) 



Moreover, if r(T) is given by then 

i(T)<^^^_, vr>i, 

2 y/i + cos2(^) 

With this to hand, we are now able to distinguish between the following cases: 

1. n < icL'r(T). In this case, the argument of the square root in ([51]) is neg- 
ative and hence both roots are imaginary. The function B{uj,n, p,T) is 
either monotonically increasing or decreasing as a function of p on [1, oo). 
Reasoning as before, note that 

f |,.,4-(T)-n. (35) 

The partial derivative vanishes precisely when n = ^ujr{T) and it is pos- 
itive for smaller n. Thus, the bound is increasing and the best possible 
bound we can find of the form (|2T]) is 



|a„| <B(w,n,l,r) = 2, n<^ujr{T). 

2. ^ajr(T) < n < . ^ u). In this case the argument of the square 

root in pip is positive. Moreover, the condition n < ujT ensures that the 
overall expression for both roots p* is positive (note that the denominator 
switches sign at n = ujT). Thus, there is a minimum of B{Lu,n, p,T) at 
p = p*(n) > 1, and this minimum satisfies p*{n) < E(T). Therefore, we 
obtain the bound 

Hence, we deduce exponential decay of the cocfEcients a„ once n exceeds 
n > ^ijjr{T). 
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3. . uj < n < wT. The arguments of the previous case still hold, 

but now p*{n) > E(T). Thus, the minimum value of B{uj,n, p,T) is 
obtained at p = E{T) and therefore 

\an\ < exp [c^Tcosh-i (2 + c(T))] . (36) 

Note that this bound is valid for all n e N, not just n in the stated range. 
However, it only gives an accurate portrait of the coefficient decay once n 
exceeds . ^ — =uj. 

^l+cos2(^) 

4. n > LuT. The denominator of ([5T|) vanishes at n = utT and switches 
sign for larger n, so that p*{n) becomes negative. Based on ([55]) . we 
conclude that the bound is monotonically decreasing as a function of p on 
[1, oo). Once more we are limited to choosing p < E{T), and therefore the 
estimate ((36)) is also applicable in this case. 



This derivation establishes the resolution constant r{T) ~ 2rsin (^) , but only 
for even oscillations f{x) ~ cosluttx. To prove the complete result, we need to 
obtain an equivalent statement for the odd functions 

f{x) = sinwTra;. (37) 

We proceed along similar lines to the cosine case. First, note that the continuous 
Fourier extension of p7p equates to the orthogonal polynomial expansion of 

sin [^rcos-i(f(i-c(r)) + i(i + c(r))] , . 
j2(uj — , == — , ue[-L,i\. 



v/i-[t(i-c(T)) + i(i + c(r))]^ 



It is not immediately apparent that this function is analytic at w = 1. However, 
simple analysis reveals that the square-root type singularity is removable, and 
consequently /2(m) is analytic in any Bernstein ellipse e{p) with p < E{T). 

Proceeding as before, we find that the nth Chebyshev polynomial coefficient 
of f2{u) admits the bound |a„| < B{uj,n, p,T), \fp > 1, where 



2 sinh 

B{uj,n,p,T) = 



.Tcosh-i [lip + i)(l - c(T)) + 1(1 + c{T))) 



p-y/Ci^ 

— 1. Suppose first that 



and Cip,T) - ^i(p+i)(l_c(r)) + i(l + c(T)) 
n < ^ujr{T). Letting p — > 1+, we find that 

|a„| < lim B{u!,n, p,T) = 2ujT. 

Next suppose that ^ujr{T) < n < . ^ to. Then, for sufficiently large w, 
we can replace B{uj, n, p, T) by 

B{uj,n,p,T) 
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T 



Figure 4: The resolution constant of Fourier extensions as a function of T. 



where B{uj, n, p, T) is the bound ([50]) of the cosine case. Hence, we deduce that 
2 



\an+k 



< 



{p*{n)Y^C{p*{n),T) 



2 n2-a;2T2 ^„ 



{p*{n)Y \j 2(1 + c{T))n^uj^T^ + {c{T)^ - l^^T* ' 

Note that the square root term is bounded since n > ^ujr{T). Indeed, it attains 
its maximal value at ti = ^u!r{T) + 1, and hence is bounded in both n and to. 
This confirms exponential decay of the coefficients a,i for n > ^ujr{T). 

In the third scenario, i.e. n > . a;, we use the value p = E{T) to 

give 

2 exp [wTcosh"^(2 + c(r))] 
Kl < -j===== . 

Note the similarity of this bound with ([55]) for the cosine case. 
To sum up, we have shown the following theorem: 

Theorem 5.2. The resolution constant r = r(T) of the continuous Fourier 
extension gn is precisely 2Tsin(^). In particular, r{T) ~ 2 for T f» 1 and 
r{T) - TT /or T > 1. 

The resolution constant r(T) is illustrated as a function of T in Fig. H] In 
Fig. [5] we confirm this theorem with several numerical examples. As discussed 
in 521 the result of the numerical computation of ([3]) may not coincide with 
the exact solution gn ■ Wc shall discuss the question of resolution power of the 
numerical solution, as opposed to the exact solution, in detail in For the 

moment, so that wc can illustrate Theorem 15.21 we have removed this issue by 
carrying out the computations in Fig. [5] in additional precision. 

As shown in the right panel of Fig. [5j the resolution constant r(T) ^ tt for 
T ^ 1. Incidentally, the fact that r(T) is independent of T for large T can be 
seen by the studying the quantity /2(u). For fixed uj, we find that 



/2(u) ==cos VT^j +C'(T-2), T^oo. 

The leading term is independent of T, thus we expect r(T) to approach a con- 
stant value as r — > oo. It is also natural to expect that this constant is the same 
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100 



(a) T = ^2 (b) T = 5 



Figure 5: The error ||/ — gn\\L^_^ , where f{x) = exp(i7ra;a;) and tu — 10, 20, 40 
(squares, circles and crosses respectively). 

as the resolution constant of polynomial approximations. Indeed, for large T 
the functions cos ^kx and sin ^kx are not at all oscillatory on [—1, 1]. In partic- 
ular, they are well approximated on [—1, 1] by their Taylor series. As such, the 
span of the first n such functions is very close to the span of their Taylor series, 
which is exactly the space of polynomials of degree n. Thus, for fixed n and 
sufficiently large T, the Fourier extension gn of an arbitrary function / closely 
resembles the best polynomial approximation to / in the sense; in other 
words, the Legendrc polynomial expansion. The value of tt for the resolution 
constant arises from the discussion in 21 

As commented at the start of this section, oscillations at 'periodic' frequen- 
cies Ld = TO € Z, are informative in that they illustrate when the Fourier 
extension behaves roughly like a classical Fourier series on [— T, T], that is, when 
T « 1, and when it does not, i.e. T ^ 1. The decay of the coefficients a„ for 
these oscillations is also quite special. Since f{x) = exp(i^7ra::) is entire and 
periodic on [— T, T], the corresponding function /2(w) given by ((26)) is also entire 
in the variable u. This means that f2{u) is analytic in any Bernstein ellipse e(r), 
and hence the corresponding Chebyshev coefficients a„ (respectively, the errors 
11/ — 5„||) decay superexponentially fast, as opposed to merely exponentially fast 
at rate E{T). Comparison of the left and right panels of Fig. [5] illustrates this 
difference. 

We now summarize this observation, along with the other convergence char- 
acteristics derived above, in the following theorem; 

Theorem 5.3. The error \\f — in approximating f{x) = cxp{iTTLUx) by 
its continuous Fourier extension gn is 0{1) for n < ^ujr{T). Once n exceeds 
^ujr{T) the error begins to decay exponentially, and when n > . ^ to, 

the rate of exponential convergence is precisely E{T), unless uj ^ ^ Jot some 
TO € Z, in which case it is superexponential. 

5.2 Points per wavelength and resolution power of the dis- 
crete Fourier extension 

The discrete Fourier extension possesses exactly the same resolution power as 
its continuous counterpart. This is a consequence of the fact that the error 



26 



committed by the polynomial interpolant of a function at Chebyshev nodes 
can be bounded by a logarithmically growing factor (the Lebesgue constant) 
multiplied by the error of its expansion in Chebyshev polynomials |35j . 

However, the discrete Fourier extension does conveniently connect the con- 
cept of resolution power with the notion of points per wavelength (see 21) • As 
shown by Thcorcm l5.2[ Fourier extensions require r{T) = 2rsin (^) points per 
wavelength to resolve oscillatory behaviour, provided the corresponding nodes 
are of the form ^T7} . 

As it transpires, the exact value for the resolution constant can also be 
heuristically obtained by looking at the maximal spacing between the mapped 
symmetric Chebyshev nodes Xi (as defined in (|17|)). Intuitively speaking, if the 
maximal node spacing for a set of n nodes is h, then one can only expect to 
be able to resolve oscillations of frequency a; < ^ . For the mapped symmetric 
Chebyshev nodes, one can show that the maximal spacing 

h = Xn ~ Xn-i —r{T), n^oo. 
Zn 

Given that there are a total of 2n + 2 nodes, one also obtains the stipulated 
value for the resolution constant in this manner. Note that identical arguments 
based on either equispaced or Chebyshev nodes in [—1,1] give maximal grid 
spacings h = ^ and h = ^ respectively. These correspond to the resolution 
constants of Fourier series and polynomial approximations, in agreement with 
previous discussions. 

5.3 Resolution power of numerical approximations 

As mentioned, it is not necessarily the case that the numerical solution of (jl4p 
(or its discrete counterpart) will coincide with the continuous (discrete) Fourier 
extension. Therefore, the predictions of Theorem 15.21 mav not be witnessed in 
computations. In Fig. |5] we give results for the numerical computation of the 
continuous/discrete Fourier extensions. For small T there is little difference with 
the left panel of Fig. [5] (except, as mentioned in ^5.3[ the computation of the 
continuous extension only attains around eight digits of accuracy). However, 
when T 3> 1 the numerically computed extension appears to possess a much 
larger resolution constant than the value r{T) ~ tt predicted by Theorem 15.21 
In fact, as demonstrated in Fig. [3 the numerical Fourier extension appears to 
have a resolution constant of 2T for all T, much like the naive estimate given 
at the start of This observation is further corroborated in Fig. [51 where a 
linear dependence of the numerical resolution constant with T is witnessed. 

This difference can be explained intuitively. The improvement of Theo- 
rem 15.21 over the nai've estimate given at the beginning of ^35. II revolves around 
the observation that slowly oscillatory functions in the frame can be recombincd 
to approximate functions with higher frequencies with exponential accuracy. 
This is due to the redundancy of the frame. However, such combinations nec- 
essarily yield large coefficients: the orthogonal polynomials that give the exact 
solution grow rapidly outside [—1,1] (see [IS]) and, hence their coefficients when 
represented in the frame consisting of exponentials that are bounded on [— T, T] 
must be large. In fact, one has the following result: 
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(a) continuous, T = \/2 



50 100 150 200 250 
(b) continuous, T = 5 




20 40 60 80 100 

(c) discrete, T = \/2 



50 100 150 200 250 

(d) discrete, T = 5 



Figure 6: The error \\f — gn\\L°° ^ ^ , where f{x) = exp(i7rajx) and u = 10,20,40. 



Lemma 5.4. Let the continuous Fourier extension g„ of f{x) — exp^iirajx) 
have coefficient vector a G C^"^"'^. Then 

\\a\\r- < ciT)E{Tr, n < K{T)uj, 
for some constant c{T) > depending on T only. 

Proof. By Parseval's relation, \\a\\i2 ~ \\gn\\L^ Writing Qn as in Theorem 
12.21 we obtain 



where a = (ao, . . . , On)^, h = (69, • ■ • 7 &n-i)^ and au and hk are given by (|9]) 
and (|10p respectively. It can be shown that 

(see [25l Thm. 3.16] for the case T = 2 - the extension to general T is straight- 
forward), and therefore 

/ n n— 1 \ 

\\ay<c[Y^\au\E{Tf + Y^ME{Tf], 



\k=0 



fe=0 



for some c > independent of n and to. By the analysis of the previous section. 



\o-k\, \bk\ ^ 1 for k < ^r{T)uj, and hence we deduce the result. 



□ 
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50 100 150 200 250 

(a) continuous 



50 100 150 200 250 
(b) discrete 



Figure 7: The error ||/ — gnWh-^-^ ^j, where f{x) = exp(50i7rx) and T = 2,3,4,5 
(squares, circles, crosses and diamonds respectively). 




5 10 20 50 100 



Figure 8: The error ||/ — ^leolkpj against cj = 1,2,..., 100, where g„ is the 
discrete Fourier extension of ,f{x) = exp(ia;7rx) and T = 2,4,8,16 (squares, 
circles, crosses and diamonds respectively). 



This lemma indicates that the coefficient vector of the continuous Fourier 
extension may well be exponentially large in the unresolved regime. As dis- 
cussed (see ^ , the numerical solution favours representations in the frame with 
bounded coefficients, since for large n the system becomes undcrdctcrmined and 
most least squares solvers will seek a solution vector with minimal norm. We 
therefore conclude that, whilst the resolution power of the frame is bounded 
by TT, the resolution power of all representations in the frame with 'reasonably 
small' coefficients is only 2T. 

In Fig.[S]we illustrate this discrepancy by comparing the 'theoretical' Fourier 
extension and its counterpart computed in standard precision. In this example, 
T = 4, i.e. E{T) « 25, and we therefore witness rapid exponential growth of 
the theoretical coefficient vector, in agreement with Lemma 15.41 In particular, 
the coefficient vector is of magnitude 10^^ at the point at which the function 
/ begins to be resolved. On the other hand, since for large n the numerical 
solver favours small norm coefficient vectors, the computed coefficient vector 
initially grows exponentially and then levels off around ||aj| w 10^^. Thus, in 
finite precision one cannot obtain the coefficient vector required to resolve / at 
the point n = ^r[T)ijj. As commented previously, these arguments can be made 
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50 100 150 200 50 100 150 

(a) Numerical extension (b) Theoretical extension 




(c) Numerical extension (d) Theoretical extension 

Figure 9: The top row shows the error |j/ — gnWh^ ^ where f{x) — cxp{iujTTx), 

(jj = 20-\/2, and T = 4. The bottom row gives the norm of the coefficient vector. 
The numerical extension was computed in Matlab, the theoretical extension was 
computed in Mathematica using higher precision arithmetic. The dotted lines 
indicate the values n = ^r{T)uj and n = ujT. 



rigorous by performing an analysis of the numerical Fourier extension [2] . 
5.4 Fixed T versus varying T 

Whilst the principal issue highlighted in the previous section, namely, that for 
large T we witness r{T) sa 2T rather than r{T) « tt in practice, is unfortunate, 
it is mainly the case T « 1 that is of interest, since this gives the highest 
resolution power. 

However, with T fixed independently of n, the resolution constant r{T) > 2; 
in other words, greater than the optimal Fourier resolution constant r = 2 
(although this difference can be made arbitrarily small). One way to formally 
obtain the value of 2 is to allow T — ^ 1+ with increasing n. For example, if 

T = l + ^, (38) 

for some c > and < a < 1, then we find that 

r(T) = 2 + — +0(7i-2")^l, n^oo, 

?7" 
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thereby giving formally optimal resolution power. The disadvantage of such a 
choice of T is that one forfeits exponential convergence. In fact, 

e(i + — )=! + — + 0{n-^'^), 

and therefore 

1 + — j - cxp(-c7rni""), n -> cx), (39) 

which indicates subexponential decay of the error (yet still spectral). Larger 
values of a lead to slow, algebraic convergence, and are consequently inadvisable. 

An alternative way to choose T is found in the literature on the KoslofF- 
Tal-Ezer mapping [51 [37] (recall the discussion in ijl.3p . In [37] the authors 
suggest choosing the mapping parameter based on some tolerance etoi. This 
limits the best achievable accuracy to 0{eto\) but gives a formally optimal time- 
step restriction. We may do the same with the Fourier extension, by solving the 
equation 

E{T)-'' = etoi, 

This gives 

T = T{n, etoi) = J (arctan (etoi)^) ■ (40) 

Note that 

rrt \ 1 log('^tol) , _2^ 

TTn 

and therefore 

frrr w o 21og(etoi) _2n 

r(T(n, etoi)) ^ 2 ^ 0(n ), n oo. 

irn 

Hence we expect formally optimal result power with this approach, as well as a 
best achievable accuracy on the order of etoi- 

As we show in the next section, choosing T in this manner works particularly 
well in practice. 



5.5 Numerical experiments 

In this final section we present numerical results for the Fourier extension applied 
to the oscillatory functions 

/i(.t) = (1 + a:2)cosl0xcosl0G7ra;, /2(a;) = Ai(-36.T - 32), (41) 

where Ai(2;) is the Airy function. Graphs of these functions are given in the top 
row of Fig. [TUl 

In the bottom row of Fig.[TO|we present the error committed by the discrete 
Fourier extension for various choices of T. The convergence rates predicted in 
the previous section are confirmed for these examples. For T decreasing with n 
the oscillations are resolved slightly sooner, but there is slower convergence in 
the resolved regime. Whether the approximation error is better or worse than 
that obtained from a fixed value of T (in this case T = |) depends on what 
level of accuracy is desired. Yet, even when high accuracy is required, a = 1/2 
(i.e. T = 1 H — Y72) appears to be a good choice for these two examples, in 
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(c) Fourier extension of /i (d) Fourier extension of /2 



Figure 10: The top row shows the functions /i and /2 defined by ((1T|) . The 
bottom row gives the error ||/ — 5„||Lp^ where T = I + ^ (squares), T = 

1 + n^i (circles), T = 1 + (crosses), T = | (diamonds) and a formaUy 
optimal T given by (^0]) with etoi =le-13 (pluses). 



spite of the theoretical subexponential convergence rate. The varying value of 
T = T{n; etoi), with tolerance set to etoi =le-13, also yields comparable results. 
The convergence rate for the case of larger a = 2/3 is slower, as expected from 
the estimate ((39| . 

As discussed, one motivation for using Fourier extensions is that, as rigor- 
ously proved in this paper, they offer a higher resolution power for small T than 
polynomial approximations, for which the resolution constant is tt (see In 
Fig. [TT] we compare the behaviour of the Chebyshev expansion and the Fourier 
extension approximation for the functions (|4T|) . Note that the latter resolves 
the oscillatory behaviour using fewer degrees of freedom, in agreement with the 
result of fj5j Indeed, for the first example function /i shown in the left panel, the 
Chebyshev expansion only begins to converge once n exceeds 150 (here, for pur- 
poses of comparison, the number of expansion coefficients is 2n), in agreement 
with a resolution constant of tt. 

The results for the second function /2 are shown in the right panel of Fig.ITT] 
The differences in resolution are smaller in this case. Yet, the fixed choice 
T = 8/7 and varying choices of T still require significantly fewer degrees of 
freedom than the Chebyshev expansion to resolve the oscillatory behaviour. As 
before, a slower asymptotic convergence rate is observed for the case. Note that 
the oscillations of /2 arc not harmonic, with the frequency increasing towards 
the right endpoint of the approximation interval. This seems to have the effect 
of reducing the advantage of schemes with optimized resolution properties. 
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50 100 150 200 50 100 150 200 

(a) Approximation of /i (b) Approximation of /2 



Figure 11: The error ||/ — gnlUpj , where g„ is Chebyshev expansion (squares), 
or the Fourier extension approximation with T = 1 + n~^/^ (circles), T — 
1 + iiT^I"^ (crosses), T = | (diamonds) and a formaUy optimal T given by (HOI) 
with etoi =le-13 (pluses). 



6 Conclusions and challenges 

The purpose of this paper was to describe and analyze the observation that 
Fourier extensions are eminently well suited for resolving oscillations. To do 
this, we derived an exact expression for the so-called resolution constant of the 
continuous/discrete Fourier extension, and offered an explanation as to why this 
constant need not be realized in the case of large T . 

Although we have touched upon the subject of numerical behaviour of Fourier 
extensions in this paper, we have not provided a full description. As mentioned, 
the continuous/discrete Fourier extensions both result in severe ill-conditioning. 
However, numerical examples in this paper and elsewhere (see [71 fTTl [30l l3T] ) 
point towards an apparent contradiction. Namely, despite this ill-conditioning, 
the best achievable accuracy can still be very high (close to machine epsilon 
for the discrete extension). It transpires that this effect can be quite compre- 
hensively explained, with the main conclusion being the following: the discrete 
Fourier extension, although the result of solving an ill-conditioned linear sys- 
tem, is in fact numerically stable. A paper on this topic is in preparation [2]. 
Incidentally, the analysis presented therein also allows one to make rigorous the 
arguments given in ij5.3l on the differing resolution constants of theoretical and 
numerical Fourier extensions for large T. 

The discrete Fourier extension introduced in this paper uses values of / on 
a particular nonequispaced mesh to compute the Fourier extension. As men- 
tioned, one important use of Fourier extensions is to solve PDE's in complex 
geometries. In this setting, such meshes are usually infeasible. For this reason, 
it is of interest to consider the question of how to compute rapidly convergent, 
numerically stable Fourier extensions with good resolution power from function 
values on equispaced meshes (as discussed, related methods based on the FE- 
Gram extension are developed in [31 [TUl 131] ) • It transpires that this can be done, 
and we will report the details in the upcoming paper [5] (see also [50]). 

Finally, we remark that there has been little discussion in this paper on the 
computational cost of computing Fourier extensions. Instead, we have focused 
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on approximation-theoretic properties, such as resolution. However, M. Lyon 
has recently introduced a fast algorithm for this purpose |31) . which allows for 
computation of Fourier extensions in O(n(logn)^) operations. 

Acknowledgements 

The authors would like to thank Rodrigo Platte for pointing out the relation 
to the Kosloff-Tal-Ezer mapping. They would also like to thank Jean-Frangois 
Maitrc for useful discussions and comments. 



References 

[1] R. A. Adams. Sobolev Spaces. Academic Press, 1975. 

[2] B. Adcock, D. Huybrcchs, and J. Martin- Vaqucro. On the numerical sta- 
bility of Fourier extensions. In preparation, 2012. 

[3] N. Albin and O. P. Bruno. A spectral FC solver for the compressible 
Navier-Stokes equations in general domains I: Explicit time-stepping. J. 
Comput. Phys., 230(16):6248-6270, 2011. 

[4] L. Badea and P. Daripa. On a fourier method of embedding domains using 
an optimal distributed control. Numer. Algorithms, 32(2-4):261-273, 2003. 

[5] F. P. Bcrtolotti, T. Herbert, and P. R. Spalart. Linear and nonlinear 
stability of the blasius boundary layer. J. Fluid Mech., 242:441-474, 1992. 

[6] J. P. Boyd. Chebyshev and Fourier Spectral Methods. Springer- Verlag, 
1989. 

[7] J. P. Boyd. A comparison of numerical algorithms for Fourier Extension 
of the first, second, and third kinds. J. Comput. Fhys., 178:118-160, 2002. 

[8] J. P. Boyd. Fourier embedded domain methods: extending a function 
defined on an irregular region to a rectangle so that the extension is spatially 
periodic and C°° . Appl. Math. Comput., 161(2):591-597, 2005. 

[9] J. P. Boyd and J. R. Ong. Exponentially-convergent strategies for defeating 
the Runge phenomenon for the approximation of non-periodic functions. I. 
Single-interval schemes. Commun. Comput. Fhys., 5(2-4):484-497, 2009. 

[10] O. Bruno and M. Lyon. High-order unconditionally stable FC-AD 
solvers for general smooth domains I. Basic elements. J. Comput. Fhys., 
229(6):2009-2033, 2010. 

[11] O. P. Bruno, Y. Han, and M. M. Pohlman. Accurate, high-order repre- 
sentation of complex three-dimensional surfaces via Fourier continuation 
analysis. J. Comput. Fhys., 227(2):1094-1125, 2007. 

[12] A. Bueno-Orovio. Fourier embedded domain methods: period and c°° 
extension of a function defined on an irregular region to a rectangle via 
convolution with Gaussian kernels. Appl. Math. Comput., 183(2):813-818, 
2006. 



34 



A. Bueno-Orovio and V. M. Perez-Garcia. Spectral smoothed boundary 
methods: The role of external boundary conditions. Numer. Methods Par- 
tial Differential Equations, 22(2):435-448, 2006. 

A. Bueno-Orovio, V. M. Perez-Garcia, and F. H. Fenton. Spectral meth- 
ods for partial differential equations on irregular domains: The spectral 
smoothed boundary method. SIAM J. Sci. Comput., 28(3):886-900, 2006. 

M. Buffat and L. Le Penven. A spectral fictitious domain method with 
internal forcing for solving elliptic PDEs. J. Comput. Phys., 230(7):2433- 
2450, 2011. 

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral meth- 
ods: Fundamentals in Single Domains. Springer, 2006. 

O. Christensen. An Introduction to Frames and Riesz Bases. Birkhauser, 
2003. 

K. S. Eckhoff. On a high order numerical method for functions with singu- 
larities. Math. Comp., 67(223):1063-1087, 1998. 

M. Garbey. On some applications of the superposition principle with 
Fourier basis. SIAM J. Sci. Comput, 22(3):1087-1116, 2000. 

D. Gottlieb and S. A. Orszag. Numerical Analysis of Spectral Methods: 
Theory and Applications. Society for Industrial and Applied Mathematics, 
1st edition, 1977. 

D. Gottlieb and G.-W. Shu. Resolution properties of the Fourier method 
for discontinuous waves. Comput. Methods Appl. Meek. Engrg, 116:27-37, 
1994. 

D. Gottlieb and G.-W. Shu. On the Gibbs' phenomenon and its resolution. 
SIAM Rev., 39(4):644~668, 1997. 

D. Gottlieb, C.-W. Shu, A. Solomonoff, and H. Vandeven. On the Gibbs 
phenomenon I: Recovering exponential accuracy from the Fourier partial 
sum of a nonperiodic analytic function. J. Comput. Appl. Math., 43(1- 
2):91-98, 1992. 

N. Hale and L. N. Trefethen. New quadrature formulas from conformal 
maps. SIAM J. Numer. Anal, 46(2):930-948, 2008. 

D. Huybrechs. On the Fourier extension of non-periodic functions. SIAM 
J. Numer. Anal, 47(6):4326-4355, 2010. 

H. H. Juang, S. Hong, and M. Kanamitsu. The NGEP regional model: An 
update. Mon. Weather Rev., 129:2904-2922, 2001. 

D. Kosloff and H. Tal-Ezer. A modified Ghebyshcv pseudospectral method 
with an 0{N~^) time step restriction. J. Comput. Phys., 104:457-469, 
1993. 

C. L. Lawson and R. J. Hanson. Solving least squares problems. Glassies 
in Applied Mathematics. SIAM, Philadelphia, 1996. 



35 



S. H. Lui. Spectral domain embedding for elliptic pdes in complex domains. 
J. Comput. Appl. Math., 225(2):541-557, 2009. 

M. Lyon. Approximation error in regularized SVD-bascd Fourier continu- 
ations. Preprint, 2011. 

M. Lyon. A fast algorithm for Fourier continuation. SI AM J. Sci. Comput., 
33(6):3241"3260, 2012. 

M. Lyon and O. Bruno. High-order unconditionally stable FC-AD solvers 
for general smooth domains IL Elliptic, parabolic and hyperbolic PDEs; 
theoretical considerations. J. Comput. Phys., 229(9):3358-3381, 2010. 

S. A. Orszag and M. Israeli. Numerical simulation of incompressible flow. 
Ann. Revs. Fluid Mech., 6:281-318, 1974. REVIEW. 

R. Pasquetti and M. Elghaoui. A spectral embedding method applied to 
the advection-difFusion equation. J. Comput. Phys., 125:464-476, 1996. 

T. Rivlin. Chebyshev Polynomials: from Approximation Theory to Algebra 
and Number Theory. Wiley New York, 1990. 

F. Sabetghadam, S. Sharafatmandjoor, and F. Norouzi. Fourier spectral 
embedded boundary solution of the Poisson's and Laplace equations with 
Dirichlet boundary conditions. J. Comput. Phys., 228(l):55-74, 2009. 

L. N. Trefethen. Is Gauss quadrature better than Clenshaw-Curtis? SIAM 
Rev., 50(l):67-87, 2008. 



36 



